Rで学ぶ実証分析
N.Yamakawa
6.1.擬似相関_1
ゲーム時間とテスト成績の散布図
data <- read.csv("data/Ch6.data.csv")
plot(data[,c(2,1)], xlab="ゲーム時間", ylab="テスト成績")
lm(grade~hours, data) %>% abline(col=4)
相関係数は-0.49です
6.1.擬似相関_2
図6.2.を再現するには「家庭環境」データが必要なため、データ生成からやり直す
plot(data[,c(2,1)], xlab="ゲーム時間", ylab="テスト成績")
lm(grade~hours, data = data0) %>% abline(col=4)
相関係数は-0.45。 テキストとちょっと違うが、許容範囲かと。
6.1.擬似相関_3
3変数間の散布図。
library(psych) #散布図用
data0 %>% pairs.panels(rug=F,ellipses=F,lm=T)
6.1.P119_図6.2的な無向グラフ
library(qgraph) #無向グラフ用
data0 %>% cor() %>% qgraph(edge.labels=T)
6.1.3.同時性
pc <- read.csv("data/police_crime.csv", stringsAsFactors=FALSE)
plot(pc[,2],pc[,3],type="n", xlim = c(1.5, 3.5),
xlab="警察官数", ylab="刑法犯認知件数")
text(pc[,2],pc[,3], pc[,1])
lm(crime~police, pc) %>% abline(col=4)
相関係数は0.129。
6.2.2.平均トリートメント効果_1
set.seed(1234); n <- 400; T <- rbinom(n, 1, 0.6)
TE <- 2; Y <- TE*T + rnorm(n); id <- 1:400
data1 <- data.frame(id, T, Y)
data1 %>% tidyr::spread(T, Y) %>% head() %>% kable()| id | 0 | 1 |
|---|---|---|
| 1 | NA | 2.485227 |
| 2 | 0.6967688 | NA |
| 3 | 0.1855139 | NA |
| 4 | 0.7007335 | NA |
| 5 | 0.3116810 | NA |
| 6 | 0.7604624 | NA |
NAの部分の値がわかればトリートメント効果が計算できるが実際はわからない。
6.2.2.平均トリートメント効果_2
boxplot(Y ~ T, data=data1)
EY1 <- mean(Y[T==1]); EY0 <- mean(Y[T==0])平均トリートメント効果(ATE)は1.872で実際の値(2)に近い。
6.2.2.平均トリートメント効果_3
set.seed(12345); n <- 400; Z <- runif(n)
T <- rbinom(n, 1, Z); TE <- 2; Y <- TE*T + (2*Z - 1) + rnorm(n)
data2 <- data.frame(T, Y, Z)
boxplot(Y ~ T, data=data2)
EY1 <- mean(Y[T==1]); EY0 <- mean(Y[T==0])平均トリートメント効果(ATE)は2.813で実際の値(2)から外れている。
7.1_外生変数
data("ToothGrowth")
summary(ToothGrowth)## len supp dose
## Min. : 4.20 OJ:30 Min. :0.500
## 1st Qu.:13.07 VC:30 1st Qu.:0.500
## Median :19.25 Median :1.000
## Mean :18.81 Mean :1.167
## 3rd Qu.:25.27 3rd Qu.:2.000
## Max. :33.90 Max. :2.000OJとVCそれぞれ30レコードずつ。 それぞれ、doseとlenの関係は？
ToothGrowth %>% lm(len~dose, data = .) -> lm.TG
ToothGrowth %>% filter(supp == "VC") %>% lm(len~dose, data = .) -> lm.VC
ToothGrowth %>% filter(supp == "OJ") %>% lm(len~dose, data = .) -> lm.OJ
summary(lm.TG)$coef## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 7.422500 1.2600826 5.890487 2.064211e-07
## dose 9.763571 0.9525329 10.250114 1.232698e-14summary(lm.VC)$coef## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 3.29500 1.427060 2.308943 2.854201e-02
## dose 11.71571 1.078756 10.860392 1.509369e-11summary(lm.OJ)$coef## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 11.550000 1.721951 6.707508 2.788784e-07
## dose 7.811429 1.301673 6.001070 1.824801e-06plot(ToothGrowth[,c(3, 1)], col = ToothGrowth$supp)
abline(lm.TG, col=4)
abline(lm.VC, col=2)
abline(lm.OJ)
VCはOJと比べて、doseが少ないときのlenの値が小さい。 doseが2を超えると、VCとOJの差は殆ど無い。
as.double(lm.TG$coef[1] + lm.TG$coef[2] *2.114)## [1] 28.06269as.double(lm.VC$coef[1] + lm.VC$coef[2] *2.114)## [1] 28.06202as.double(lm.OJ$coef[1] + lm.OJ$coef[2] *2.114)## [1] 28.063367.2_内生変数
set.seed(1234)
n <- 200; e <- rnorm(n)
X <- (1+0.5*e)*runif(n)
b0 <- 1; b1 <- 2
Y <- b0 + X*b1 + e
lm(Y~X)$coef## (Intercept) X
## 0.02967807 3.71966851Xの係数3.7196685はb1の推定としては不適切？
cov(X,e)/var(X) + b1## [1] 3.719669数式(7.7)による補正rをすると、上記Xの係数とほぼ一致
7.3.2_双子データ
OpenMxパッケージの双子データを用いる
library(OpenMx)
data(twinData);head(twinData)## fam age zyg part wt1 wt2 ht1 ht2 htwt1 htwt2 bmi1 bmi2
## 1 115 21 1 2 58 57 1.7000 1.7000 20.0692 19.7232 20.9943 20.8726
## 2 121 24 1 2 54 53 1.6299 1.6299 20.3244 19.9481 21.0828 20.9519
## 3 158 21 1 2 55 50 1.6499 1.6799 20.2020 17.7154 21.0405 20.1210
## 4 172 21 1 2 66 76 1.5698 1.6499 26.7759 27.9155 23.0125 23.3043
## 5 182 19 1 2 50 48 1.6099 1.6299 19.2894 18.0662 20.7169 20.2583
## 6 199 26 1 2 60 60 1.5999 1.5698 23.4375 24.3418 22.0804 22.3454
## cohort zygosity age1 age2
## 1 younger MZFF 21 21
## 2 younger MZFF 24 24
## 3 younger MZFF 21 21
## 4 younger MZFF 21 21
## 5 younger MZFF 19 19
## 6 younger MZFF 26 26fam：家族番号
age：双子の年齢
zyg：遺伝子とコホートをコード化(1-10)
part：？(0:15, 1:224, 2:3569)
wt1,2：双子の体重(kg)
ht1,2：双子の身長(m)
htwt1,2：双子のBMI(kg/m^2)
bmi1,2：双子のBMIの対数
cohort：youngerかolder
zygosity：遺伝子パターン(MZFF, MZMM, DZFF, DZMM, DZOS)
boxplot(age~cohort, twinData)
cohortの境界は30歳
双子の組に対して、身長の差と体重の差で線形モデルを作る
twinData %>%
mutate(wt0 = wt1 - wt2, ht0 = ht1 - ht2) -> data1
lm(wt0 ~ ht0-1, data=data1) -> lm.data1
summary(lm.data1)$coef## Estimate Std. Error t value Pr(>|t|)
## ht0 87.56072 1.661656 52.69486 0plot(data1$ht0, data1$wt0, xlab="身長", ylab="体重")
abline(lm.data1, col=4)
身長1cm毎に体重は0.87kg増える
c("age", "zyg", "part", "wt0", "ht0", "cohort", "zygosity") -> coln
rbind(
twinData[,c(2:4,5,7,13:14)] %>% set_colnames(coln),
twinData[,c(2:4,6,8,13:14)] %>% set_colnames(coln)
)-> data2
lm(wt0 ~ ., data=data2) -> lm.data2
summary(lm.data2)$coef## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -67.0140567 2.42647491 -27.6178651 7.880592e-160
## age 0.1270048 0.01097231 11.5750260 1.017459e-30
## zyg 0.6412055 0.07149348 8.9687268 3.752456e-19
## part 0.1472285 0.57012365 0.2582397 7.962292e-01
## ht0 72.6191300 1.31619664 55.1734654 0.000000e+00
## cohortyounger 1.0833814 0.46770506 2.3163773 2.056505e-02
## zygosityMZMM 3.0815594 0.32938609 9.3554632 1.084290e-20
## zygosityDZFF -1.0652820 0.26069244 -4.0863555 4.428643e-05
## zygosityDZMM 2.4107945 0.37732550 6.3891639 1.769711e-10身長1cm毎に体重は0.72kg増える。